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ABSTRACT 

We present a three-dimensional map of interstellar dust reddening, covering three-quarters of the 
sky out to a distance of several kiloparsecs, based on Pan-STARRS 1 and 2MASS photometry. The 
map reveals a wealth of detailed structure, from filaments to large cloud complexes. The map has 
a hybrid angular resolution, with most of the map at an angular resolution of 3.4' to 13.7', and a 
maximum distance resolution of ~25%. The three-dimensional distribution of dust is determined in a 
fully probabilistic framework, yielding the uncertainty in the reddening distribution along each line of 
sight, as well as stellar distances, reddenings and classifications for 800 million stars detected by Pan- 
STARRS 1. We demonstrate the consistency of our reddening estimates with those of two-dimensional 
emission-based maps of dust reddening. In particular, we find agreement with the Planck 7353 c;iwr 
based reddening map to within 0.05 mag in E (B — V) to a depth of 0.5 mag, and explore systematics 
at reddenings less than E(B — V) ~ 0.08 mag. We validate our per-star reddening estimates by com¬ 
parison with reddening estimates for stars with both SDSS photometry and SEGUE spectral classifi¬ 
cations, finding per-star agreement to within 0.1 mag out to a stellar E (B — V) of 1 mag. We compare 
our map to two existing three-dimensional dust maps, by Marshall et al. (2006) and Lallement et al. 
(2013), demonstrating our finer angular resolution, and better distance resolution compared to the 
former within ~3kpc. The map can be queried or downloaded at http://argonaut.skymaps.info. 

We expect the three-dimensional reddening map presented here to find a wide range of uses, among 
them correcting for reddening and extinction for objects embedded in the plane of the Galaxy, studies 
of Galactic structure, calibration of future emission-based dust maps and determining distances to 
objects of known reddening. 


1. INTRODUCTION 

The Milky Way is the only Galaxy that can be observed 
in such close detail, yet most of the plane of the Galaxy 
is veiled by dust. While dust makes up just 1% of the 
mass of the interstellar medium, it absorbs on the order 
of 30% of all starlight (Draine 2003). Interstellar dust 
absorbs and scatters light in the ultraviolet (UV), opti¬ 
cal and near-infrared (NIR), re-radiating it in the mid- to 
far-infrared (MIR and FIR, respectively). Understand¬ 
ing the spatial distribution of dust is therefore critical 
for UV, optical and NIR astronomy, where dust is an ex¬ 
tinguishing foreground, to extragalactic astronomy and 
cosmology, where it is a radiating foreground, as well as 
to star formation, where the dust itself is a primary ob¬ 
ject of study. Detailed studies of stellar populations and 
substructures in the Galaxy require accurate corrections 
for extinction and reddening due to dust. The plane of 
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the Milky Way, containing the majority of the Galaxy’s 
stellar content, is also the most heavily extinguished re¬ 
gion. But dust is not only a nuisance to astronomers. 
Dust traces the interstellar medium, and a clearer pict ure 
of the spatial distribution of dust would aid in under¬ 
standing the processes that shape the Galaxy, from star 
formation to how feedback from supernovae and stellar 
winds shape our Galaxy’s ISM. 

Dust can be mapped either through its extinction or 
through its emission. The spectral energy distribution 
and amplitude of FIR dust emission is most sensitive 
to dust column density, temperature, and grain-size dis¬ 
tribution. By modeling these properties across the sky, 
a map of dust column density may be obtained, which 
can then be converted to extinction or reddening by 
calibration against astrophysical extinction or redden¬ 
ing standards. Such emission-based methods recover 
the angular distribution of dust, but not its distribu¬ 
tion in distance. Burstein & Heiles (1982) produced an 
all-sky map of dust reddening based on III emission and 
galaxy counts, assuming the gas and dust to be uniformly 
mixed. Schlegel et al. (1998) (hereafter “SFD”) produced 
a widely-used all-sky reddening map, using DIR BE and 
IRAS maps of FIR emission to model dust column den¬ 
sity and temperature. More recently, the Planck Collab¬ 
oration has released all-sky reddening maps based on a 
similar emission-modeling technique (Planck Collabora¬ 
tion et al. 2014). 

A second class of dust maps is based on extinction 
or reddening estimates of sources distributed across the 
sky. Lada et al. (1994, NICE) compared the average 
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II — I\ colors and number counts of stars in target and 
control fields to map dust reddening in two dimensions. 
Lombardi & Alves (2001. NICER) and Lombardi (2009, 
NICEST) extended on this algorithm, which has now 
been applied to 2MASS photometry to obtain redden¬ 
ing maps of a host of cloud complexes (Lombardi et al. 
2006, 2008, 2010, 2011; Alves et al. 2014). Froebrich 
applied a similar color excess method to 2MASS pho¬ 
tometry to produce a two-dimensional, all-sky reddening 
map (Rowles & Froebrich 2009). Schlafly et al. (2010) 
used the location of main-sequence turn-off stars in color- 
color space to simultaneously map dust reddening and 
test different reddening laws. Peek & Graves (2010) used 
passively evolving galaxies as a standard color source in 
order to map reddening and correct the SFD map. 

Because stars are distributed throughout, the Galaxy, 
they can be used to trace the dust distribution in three di¬ 
mensions. Most methods relying on this principle group 
stars into separate sightlines, and then determine stellar 
reddening as a function of distance along each sighthne. 
The challenge with such methods is to simultaneously 
determine stellar type (and thus the intrisic stellar col¬ 
ors and luminosity), distance and reddening on the basis 
of photometry alone. Marshall et al. (2006) developed 
a method that iteratively improves distance and redden¬ 
ing estimates to post-main sequence stars, updating the 
dust column in each distance bin with each iteration so 
that the intrinsic stellar colors match those predicted by 
the Besangon model of the Galactic stellar population 
(Robin et al. 2003). Marshall et al. (2006) applied this 
method to 2MASS photometry of the Galactic plane, 
producing a three-dimensional reddening map out to a 
distance of several kiloparsecs in the region |5| < 10°, 

— 100° < l < 100°. A number of groups have pursued 
methods that determine maximum-likelihood parameters 
for the stars along each sightline, combining the esti¬ 
mates in each sightline into a distance-reddening rela¬ 
tion. Berry et al. (2011) applied such a method to SDSS 
photometry, producing a 3D reddening map and measur¬ 
ing variation in the dust extinction spectrum across the 
SDSS footprint. Chen et al. (2014) used optical XSTPS- 
GAC, NIR 2MASS and MIR WISE photometry to pro¬ 
duce a 3D reddening map of the Galactic anticenter. Sale 
(2012) developed a probabilistic framework to simulta¬ 
neously infer stellar parameters and the dust extinction 
distribution along each line of sight, and Sale et al. (2014) 
applied this method to IP1IAS photometry of the north¬ 
ern Galactic plane. Hanson & Bailer-Jones (2014) devel¬ 
oped a similar probabilistic framework for inferring stel¬ 
lar types, distances and dust foreground properties, and 
applied this method to SDSS and UKIRT Infrared Deep 
Sky Survey photometry in regions of the high-Galactic- 
latitude sky. Lallement et al. (2013) took a somewhat 
different approach, using ~ 23,000 stellar parallaxes and 
reddening estimates from a number of sources to infer the 
3D distribution of dust opacity out to a distance of 800 

- 1000 pc in the plane of the Galaxy, and ~ 300 pc out 
of the plane. Green et al. (2014) developed a probabilis¬ 
tic method for determining dust reddening in 3D from 
stellar photometry, and Schlafly et al. (2014c) presented 
a map of dust reddening integrated to 4.5 kpc, covering 
three quarters of the sky, based on applying this method 
to optical Pan-STARRS 1 stellar photometry. Using this 
same 3D mapping technique, Schlafly et al. (2014a) found 


that the Orion molecular complex appears to form part 
of a larger bubble structure. Schlafly et al. (2014b) deter¬ 
mined the distance to a large number of molecular clouds 
by applying a related method to Pan-STARRS 1 stellar 
photometry. 

In this paper, we present a three-dimensional map 
of dust reddening over three quarters of the sky. Our 
map covers the Northern sky down to S ~ —30°, has 
an adaptive resolution of typical angular scale 3.4' to 
13.7', logarithmically-spaced distance bins with a width 
of ~25%, and extends to E (D — V) ~ 1.5 — 2mag. Our 
map is based on high-quality Pan-STARRS 1 (PS1) pho¬ 
tometry for 800 million stars, ~ 200 million of which 
have matched 2MASS photometry. The reddening along 
each sightline is determined in a fully probabilistic man¬ 
ner, according to the method developed in Green et al. 
(2014). We strongly encourage readers who wish to un¬ 
derstand our method for inferring the three-dimensional 
distribution of dust reddening to read that paper, as the 
present paper focuses primarily on results obtained using 
that method. 

This paper is organized as follows. In §2, we briefly 
describe the method developed in Green et al. (2014), 
and describe the priors we place on the three-dimensional 
distribution of dust reddening. In §3, we describe the 
data sources which go into our map. In §4, we discuss 
results from the map, and in §5, we compare our map 
to existing 3D and 2D dust maps. We describe how the 
map can be accessed in §6. 

2. METHOD 

Our three-dimensional dust map is constructed pixel- 
by-pixel from independently determined line-of-sight red¬ 
dening profiles. We group stars into sightlines with typ¬ 
ical angular scales of ~ 6.8', probabilistically inferring 
the reddening, distance and stellar type of each star in¬ 
dependently. We then assume that the distances and 
reddenings of the stars in a given line of sight he along 
a single profile, with reddening increasing monotonically 
with distance. We sample from the posterior density of 
distance-reddening profiles, returning the uncertainty in 
the line-of-sight reddening in the pixel. By repeating this 
process in each angular pixel, we obtain the 3D distribu¬ 
tion of dust throughout the Pan-STARRS 1 survey vol¬ 
ume (PSl; Kaiser et al. 2010). The details of our method 
are described in Green et al. (2014). 

We have made three notable modifications to our 
method since publication of Green et al. (2014). The 
first modification is the inclusion of stellar photometry 
from the Two Micron All Sky Survey (2MASS). We com¬ 
pile joint PS1+2MASS stellar templates, in a manner 
very similar to the process laid out for the PSl templates 
used in Green et al. (2014). Appendix A contains details 
about the construction of our PS1+2MASS templates, 
while Appendix B details our estimate of the 2MASS 
survey depth. 

The second modification to our method is to allow the 
reddening of individual stars to deviate slightly from the 
overall line-of-sight reddening profile in a pixel. Our ba¬ 
sic model assumes that within an angular pixel, the dust 
density does not depend on angle, but is solely a function 
of distance. This is a good assumption at very fine an¬ 
gular scales, but it must obviously begin to break down 
at some angular scale, depending on the power spectrum 
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of the dust density field. Because each star in an an¬ 
gular pixel is behind a slightly different column of dust, 
we allow the total column of dust in front of each star 
to vary slightly from the modeled column density. For a 
pixel with a small angular scale, in which the dust col- 
umn should be more uniform across the pixel, the scatter 
of the stars around the modeled distance-reddening pro¬ 
file should be small. For pixels with larger angular scale, 
where angular variations in the dust column are correp- 
sondingly greater, one should allow the stars to deviate 
from the modeled distance-reddening profile to a greater 
extent. We therefore include an additional parameter for 
each star, describing the fractional variation of the star 
from the modeled distance-reddening profile in its pixel. 
This modification to the model amounts to introducing a 
smoothing factor in the probability density of each star in 
distance-reddening space, and therefore contributes only 
negligibly to the overall computation time. We discuss 
this modification in more detail in§2.2. 

A third modification we have made is to adopt the 
increased, bias-corrected stellar disk scale lengths and 
heights presented in Table 10 of .Juric et al. (2008) for 
our disk priors. We also adopt a less dense halo prior 
than Juric et al. (2008), reducing the local halo density 
(relative to the local thin-disk stellar density) from /;, = 
0.0051 to /;, = 0.0006, a similar value to that used by 
the Besangon model (Robin et al. 2003). 

A brief summary of our model, including the changes 
introduced in §2.2, is given is §2.3. 

Finally, we improve the convergence of our Markov 
Chain Monte Carlo sampling of the line-of-sight redden¬ 
ing profile by introducing a new type of proposal step, 
which we term the “swap” proposal. See Appendix C for 
details. 

2.1. 3D Dust Priors 

Green et al. (2014) describes how we use photometry of 
stars to determine the distribution of dust in the Galaxy. 
We begin by separating the sky into small sightlines. In 
a small region of the sky across which the dust column 
does not vary much angularly, more distant stars should 
be more heavily reddened than nearby stars, as all the 
stars he along the same dust column. We therefore begin 
by calculating a probabilistic reddening and distance es¬ 
timate for each star in the sightline. We then determine 
the amount of dust reddening in each distance bin along 
the sightline, under the constraint that the line-of-sight 
distance vs. reddening profile has to be consistent with 
our distance and reddening inferences for all the stars in 
the sightline. 

Green et al. (2014) states our map-making method gen¬ 
erally, without choosing specific priors on the distribu¬ 
tion of reddening in 3D. Here, we will briefly restate the 
formalism used in Green et al. (2014), and then define 
our 3D reddening priors. 

Our goal of inferring the line-of-sight dust reddening on 
the basis of stellar photometry begins with inferring the 
distance and reddening of each star based on its photom¬ 
etry. Label each star in the line of sight by a number i, 
and denote the probability density for an individual star 
to he at distance modulus, p, and reddening, E, given 
photometry m, by 

p ( m , Ei \ fhi ) . (1) 


We precompute p{pi , E; \ rhi ) using a kernel density esti¬ 
mate of Markov Chain Monte Carlo samples. This tech¬ 
nique, similar to that used in Hogg et al. (2010), is not 
guaranteed to converge to exactly the target density. It 
is possible, however, to substitute a mathematically cor¬ 
rect, but somewhat more computationally expensive al¬ 
gorithm, pseudo-marginal Markov Chain Monte Carlo, 
(Beaumont 2003; Andrieu & Roberts 2009), in which a 
new noisy estimate of a marginalized likelihood term is 
computed at each Markov chain step. However, we ex¬ 
pect any possible bias introduced by precomputing Eq. 
(l) once to be small due to the relatively large number 
of stars in each sightline. Compared to possible inaccu¬ 
racies in our stellar model, for example, we expect the 
inexactness of our sampling method to be a small effect. 

As explained in Green et al. (2014), we place a flat 
prior on E, for each individual star when precomputing 
Eq. (1), as the prior on reddening will come in when we 
combine information from all of the stars in a sightline. 
We now parameterize the line-of-sight reddening by a set 
of parameters, a, so that the cumulative reddening out 
to a distance modulus p can be written as: 

a). (2) 

The probability density of <5 is then given by a product 
over line integrals through the individual stellar proba¬ 
bility density functions, following E(p\a): 

p(a | {m}) oc p{a) ^d//., p(p u E{pp a) | m,) . (3) 

We parameterize the reddening in each line of sight 
as a monotonically increasing, piecewise-linear function 
in distance modulus. Dividing up the line of sight into 
bins of equal width in distance modulus, we sample the 
increase in reddening in each bin, so that 

a = {AE k \k = l,2,...,N hins } , (4) 

where k denotes the bin index, and A r bi ns is the number 
of distance bins. We place a log-normal prior on the 
reddening in each bin: 

A E k = e ek , where e k ~ Jf(e k , <r E ) . (5) 

In each bin, we set e k so that the mean reddening in 
each voxel matches what would be expected from the 
smooth disk component of Drimmel & Spergel (2001), 
up to a constant normalization that is the same for the 
entire map. In order to fully define the priors, there are 
therefore two global parameters that must be set: 

1. dE(B—vj/dsl the local normalization of the dust 

' ls=0’ 

reddening per unit distance, and 

2. (7 e , the width of the log-normal prior on dust col¬ 
umn density in each voxel. 

We fit these two numbers so that our priors, in projec¬ 
tion, produce two-dimensional dust maps with the same 
overall normalization and standard deviation as the SFD 
dust map. We find that <r f = 1.4 and a local reddening 
per unit distance of = 0.2magkpc _1 roughly 

matches the mean and variance of the SFD dust map 
over large angular scales. 

Fig. 1 shows a random realization of our 3D dust pri¬ 
ors, next to the SFD dust map. In each bin, we limit e 
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SFD Priors 



Fig. 1. The SFD dust map, compared to a random draw from the priors on the 3D dust reddening distribution, used in the construction 
of our map. The top two panels, from left to right, show the SFD reddening and the 2D reddening map that results from a draw from the 
priors, both on a log scale. The bottom two panels show the same maps, after dividing out the mean projected reddening in the priors. 
The priors are limited so that the mean expected reddening in any given distance bin does not exceed a pre-defined amount. This is done 
in order to avoid inferring large amounts of reddening in the absence of data. Regions that are affected by this clipping are shaded in gray. 
The priors do not include spatial correlations, as can be seen by comparing the bottom two panels. 


to the range —12 < l < —4. On the lower end, the limit 
helps our fit to converge by preventing the prior from 
becoming too stringent. On the upper end, the limit 
prevents our fit from inferring large amounts of dust in 
the absence of data. Regions where e has been limited 
to —1 are shaded in Fig. 1. Our priors do not impose 
spatial correlations across lines of sight, and thus in the 
absence of data, produce cloud-free maps. The detailed 
cloud structures that emerge from our 3D dust modeling 
therefore derive entirely from the data, rather than from 
the priors. 

It is worth noting that our priors diverge from what 
one should expect for the real distribution of dust in a 
number of ways. In order to render the problem of fit¬ 
ting the 3D distribution of dust tractable, we fit the dust 
column along each sightline separately. We know, how¬ 
ever, that in real life, dust density has spatial correla¬ 
tions. Moreover, there is nothing special about the Sun’s 
place in the Galaxy, but our dust map voxelizes the sky 
into pencil beams centered on the Solar System. This of 
course makes practical sense, since angles are much easier 
to measure than distances in astronomy, but it is again 
an unreal feature of our voxelization. A more realistic 
model would treat the dust density as a continuous field, 
or voxelize the Galaxy in a way that treats the angular 
and radial directions equally, and would impose corre¬ 
lations between dust density in nearby points in space 
(see, e.g., Lallement et al. 2013; Sale & Magorrian 2014). 
This entails significantly more algorithmic and compu¬ 
tational complexity than the method used here, and we 
defer such work to the future. Within the constraints of 
our present setup independent sightlines with pencil- 
beam-like voxels our priors attempt to reasonably trace 
the properties of the Galaxy, including the mean dust 
density in each voxel and the overall variance in dust 
column across the sky. 

2.2. Scatter in Line-of-Sight Reddening Profile 


A basic assumption of our model, as described in Green 
et al. (2014), is that within a given pixel, the dust density 
varies only with distance, and not with angle. If the 
pixels are sufficiently small, this is a good assumption. 
We are limited, however, in how small we can make each 
pixel by the need to include enough stars in each pixel to 
probe the dust density at a range of distances. Increasing 
the angular resolution of the map decreases the number 
of stars in each pixel, effectively decreasing the distance 
resolution of the map. We have found that we obtain 
best results for pixels containing a few hundred stars, 
and we vary the resolution of our pixels across the sky to 
obtain approximately the same number of stars in each 
pixel (see Fig. 2). 

Note that varying the pixel size across the sky in this 
way technically violates the principles of Bayesian infer¬ 
ence. From a forward-modeling point of view, the dis¬ 
tribution of dust influences the number of stars that are 
observed in any given region of the sky. Sale (2015), 
for example, discusses how a catalog of stars observed 
in a survey can be described as a Poisson point process 
whose rate is determined by the distribution of stars in 
the Galaxy, the three-dimensional distribution of dust, 
and the survey selection function. Yet we are using the 
number of stars observed in each part of the sky to deter¬ 
mine how to pixelize the sky. Our pixehzation is therefore 
set, in part, by an observable (the density of stars across 
the sky) that is a consequence of the model. Neverthe¬ 
less, as we are treating each pixel independently, and we 
would like to keep pixels as small as possible without re¬ 
ducing the number of stars per pixel below a few hundred, 
this rules violation is difficult to avoid. We expect the 
impact of this violation to be small in most regions of the 
sky. In regions with large sub-pixel variation in dust red¬ 
dening, however, our map may be biased towards lower 
reddenings, as we preferentially detect stars in regions of 
the pixel with lower reddening. 

The typical resolution of the map is 6.8' x 6.8', cor- 
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Fig. 2. Assigning stars to different pixels, in which the line- 
of-sight dust profile will be fit independently. Here, each dot rep¬ 
resents a point-source detection, and the dots are colored by the 
HEALPix pixel they are assigned to. The pixel scale varies in or¬ 
der to keep the number of stars per pixel roughly constant. Note 
that some pixels show significant extinction-induced variations in 
stellar density. 

responding to an nside = 512 HEALPix pixelization 
(Gorski et al. 2005). At this resolution, there can still 
be significant power in the dust density spectrum below 
the pixel scale. This can pose problems for our method, 
especially in the vicinity of dense clouds and filamen¬ 
tary structures, where the sub-pixel angular variation is 
largest . In order to deal with sub-pixel angular variation 
in the dust density, we relax our assumption that all stars 
lie along the same dust column by allowing each star to 
deviate from the local “average” dust column by a small 
amount . The reddening of star i is parameterized as 

Ei = (1 + 5,) £(//,;«) , (6) 

where S, is the fractional offset of the star from the local 
dust column. E(pj\d). 

In effect, our model is therefore that within each 
HEALPix pixel, the reddening is a white noise process, 
with a mean that increases piecewise linearly with dis¬ 
tance. Each star samples this white noise process at a 
particular distance and angular position within the pixel. 
The parameter S, is then understood as the fractional 
residual (from the mean reddening in the pixel at the 
given distance) of the reddening column at the angular 
location and distance of star i. 

We put a Gaussian prior on S, , with zero mean and 
standard deviation dependent to the scale of the pixel 
(allowing more variation in larger pixels) and the local 
dust column (allowing greater fractional variation in re¬ 
gions of greater reddening). In detail, 

p(6i \m, a) =M(5i 10, a 5 ) , (7) 

with 

as = ciE{p.i\ a) + b. (8) 

Here, a and b are parameters that we set in order to 
match the variation we see at the given pixel angular 
scale in the Planck radiance-based two-dimensional dust 


map. We compute the RMS scatter within HEALPix 
pixels of difference scales, finding that the scatter is well 
described by setting the coefficients a and b to 

logic a = 0-88 log 10 ( p) - 2-90 , (9) 

log 10 6 = 0.58log 10 (p) -1-88, (10) 

where is the angular pixel scale, defined as the square- 
root of the pixel solid angle. 

We have found through trial and error that it is prefer¬ 
able to impose a minimum scatter of ~ 10% on the in¬ 
pixel dust col umn . Additionally, if we allow 6, to ap¬ 
proach unity, we clearly risk the possibility of scattering 
a star to negative dust column. We therefore never allow 
(7,5 > 0.25. 

Although the complication introduced in this section 
adds an additional parameter, Sj, for each star, it can 
be achieved with minimal additional computational re¬ 
sources. We are introducing a Gaussian scatter in the 
reddening of each star from the “average” reddening 
in the pixel, and an appropriate Gaussian smoothing 
of the individual stellar probability density surfaces, 
p(pi, E, | m,), achieves this effect. Appendix D de¬ 
tails how the individual stellar probability surfaces are 
smoothed. 

2.3. Summary of Model 

In summary, our model of each sightline contains the 
following elements: 

• The increase in the “average” reddening in each 
distance bin: a = {A Ek \ k = 1,..., nbi ns }. 

• The distance modulus, , stellar type, ©, and frac¬ 
tional offset, Sj from the “average” reddening of 
each star i. 

The reddening of star i is therefore determined both by 
a and 6 ,. Together with the distance and type of the 
star, one can obtain model apparent magnitudes for the 
star, and thus a likelihood: 

|/q,0i,<5i,o) . (11) 

We also have per-star priors on distance and stellar type, 
given by a smooth model of the distribution of stars 
throughout the Galaxy, and a prior on the offset of each 
star from the local reddening column: 

p(p,i,Qi, h | a) = pi^ii ©i) p{&i | <*) • (12) 

We finally have a log-normal prior on the increase in 
reddening in each distance bin along each sightline, p(a), 
whose mean is chosen to match a smooth model of the 
distribution of dust throughout the Galaxy. 

The posterior on the line-of-sight reddening along one 
sightline is then given by 

"stars . 

p(a | {in}) = p{a) p[rhi \ p u ©i, Si, a) 

1=1 

x p(Si\pi,d) . (13) 
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We pre-compute the likelihood and prior terms for the 
individual stars, and then sample in a. The details of 
how this is done are given in Green et al. (2014), and in 
Appendix D of this paper. 


that they correspond to the model sketched out directly 
above. As shown in detail in Appendix E, the correct 
reweighting of the stellar Markov chain samples, assum¬ 
ing a particular line-of-sight reddening profile a, is 


2.4. Improving Stellar Inferences using the 3D Map 

In order to create the 3D dust map, we first proba¬ 
bilistically infer the distance and reddening to each star 
individually, with no information about the 3D structure 
of the dust. After creating the 3D dust map, however, we 
have a very strong constraint on how reddening should 
increase with distance, and this should impact our infer¬ 
ences about individual stellar distances and reddenings. 
In order to improve our stellar parameter inferences, we 
should replace our initial assumption about stellar red¬ 
dening (a flat prior) with a new one that favors stellar 
reddenings close to the measured distance-reddening re¬ 
lation along the sightline. This can be done, in practice, 
by reweighting the probability density functions we ini¬ 
tially calculated for each star. In this section, we there¬ 
fore define a reweighting of the Markov Chain samples 
for the individual stars which takes into account the line- 
of-sight reddening. 
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Fig. 3.- Reweighting Markov Chain samples for an individ¬ 
ual star, based on the line-of-sight reddening profile. Each black 
line is a reddening profile drawn from the posterior on line-of-sight 
reddening. The dots are posterior samples of parameters (includ¬ 
ing distance and reddening) for one star viewed in isolation, i.e., 
not conditioned on the line-of-sight reddening. The samples are 
reweighted in order to condition them on the line-of-sight redden¬ 
ing, assigning greater weight to samples that are consistent with 
the reddening profile. Including information about the line-of-sight 
reddening can significantly reduce the uncertainties in stellar dis¬ 
tance, reddening and type. 


Let. us first consider the case in which we fix the line- 
of-sight reddening profile. In the formalism used here, 
that means that we fix the parameters ci, expressing the 
line-of-sight reddening profile as E(p;a). The reddening 
of an individual star is then determined by the distance 
modulus of the star, /t, as well as the fractional deviation, 
4 of the stellar reddening from the local dust column. As 
before, there are also parameters representing the stellar 
type, 0. In this formulation, the posterior density of 
the stellar parameters, given its photometry in and the 
line-of-sight reddening profile, is determined by 

p(p, S | m, a) oc p(rh \ p, 0, <5, cv) p[p, ©) p(5 | //., a) . 

(14) 

We already have Markov chain samples in distance, red¬ 
dening and stellar type for each star, for a model which 
does not take the line-of-sight reddening into account. 
We would like to apply weights to these samples, so 


lL’k <x 


A/~(4-|o^,s) 

E{g.k-,a) 


(15) 


where k indexes the sample, and 4 is the fractional off¬ 
set of the sample reddening, £*, from the line-of-sight 
reddening, E(pk', <*), at distance pk- The numerator in 
the above weight corresponds to a prior on the offset of 
the stellar reddening from the local dust column, while 
the denominator comes from the Jacobian transforma¬ 
tion from reddening to fractional offset from the local 
reddening. 

Reweighting each of the stellar parameter samples by 
the above factor, and normalizing the sum of the weights 
to unity, we obtain an inference for the stellar parame¬ 
ters, conditioned on a particular line-of-sight reddening 
profile, E(p:a). We marginalize over the line-of-sight 
reddening profile by repeating this procedure for each 
sampled value of a, summing the weight applied to each 
stellar sample. 

It. might be objected that since we infer a using the 
photometry of all the stars in the given pixel, the sam¬ 
ples we drew for a in our initial processing are already 
dependent on the photometry for the star whose sam¬ 
ples we are reweighting. Put more formally, the prior 
we use here when marginalizing over line-of-sight red¬ 
denings is the inference we obtained earlier, p(a | {?n}), 
which is conditional on all the photometry in the pixel, 
{m}. We are using the photometry of a given star to 
infer the line-of-sight reddening, and then again to infer 
the stellar parameters, conditional on that line-of-sight 
reddening. 

This would be a problem if the photometry of any given 
star significantly affected the line-of-sight reddening in¬ 
ference. However, if photometry from a large number of 
stars informs the line-of-sight reddening, then no one star 
should have a significant impact on the inferred line-of- 
sight reddening profile. In a more formally correct for¬ 
mulation of the problem, we would first infer the line-of- 
sight reddening using all the stars in the pixel except the 
star whose parameters we wish to infer, and we would 
then infer the parameters for that star, conditional on 
the inferred line-of-sight reddening profile. As each pixel 
contains hundreds of stars, however, we expect our pro¬ 
cedure to approximate this formally correct procedure 
closely. 

Fig. 3 shows the result of reweighting the samples for 
one star. Depending on the line-of-sight reddening profile 
and the distribution of the unweighted stellar parameter 
samples, individual stellar inferences can be tightened 
dramatically by taking the line-of-sight reddening into 
account. Our knowledge of the parameters describing 
one star is therefore dependent not only on its photom¬ 
etry, but also on the photometry of its neighbors on the 
sky. Because neighboring stars lie along the same dust 
column, inferences for nearby stars are coupled through 
the requirement that the dust column increase with dis¬ 
tance. 


3. DATA 











3.1. Pan-STARRS 1 

Pan-STARRS 1 is a 1.8-meter optical and near-infrared 
telescope located on Mount Haleakala, Hawaii (Kaiser 
et al. 2010; Ilodapp et al. 2004). The telescope is 
equipped with the GigaPixel Camera 1 (GPC1), consist¬ 
ing of an array of 60 CCD detectors, each 4800 pixels on 
a side (Tonry k. Onaka 2009; Onaka et al. 2008). From 
May 2010 to April 2014, the majority of the observing 
time was dedicated to a multi-epoch 3 tt steradian sur¬ 
vey of the sky north of 6 = —30° (Chambers in prep.). 
The 37 t survey observes in five passbands gpi, rpi, i p i. 
zpi, and ypi, similar to the Sloan Digital Sky Survey 
(SDSS; York et al. 2000), with the most significant dif¬ 
ference being the replacement of the Sloan u band with 
a near-infrared band, yp\. The PS1 filter set spans 400 
1000 nm (Stubbs et al. 2010). The images are processed 
by the Pan-STARRS 1 Image Processing Pipeline (IPP) 
(Magnier 2006), which performs automatic astrometry 
(Magnier et al. 2008) and photometry (Magnier 2007). 
The data is photometrically calibrated to better than 1% 
accuracy (Schlafly et al. 2012; Tonry et al. 2012). The 3 tt 
survey reaches typical single-exposure depths of 22 mags 
(AB) in gpi, 21.5 mags in rpi and ipj, 20.8 mags in zpi, 
and 20 mags in yp\. The resulting homogeneous optical 
and near-infrared coverage of three quarters of the sky 
makes the Pan-STARRSl data ideal for studies of the 
distribution of the Galaxy’s dust. 

3.2. Two Micron All Sky Survey 

The Two Micron All Sky Survey (2MASS) is a uni¬ 
form all-sky survey in three near-infrared bandpasses, 
.7, H and K„ (Skrutskie et al. 2006). The survey de¬ 
rives its name from the wavelength range covered by the 
longest-wavelength band, K s , which lies in the longest- 
wavelength atmospheric window not severely affected by 
background thermal emission (Skrutskie et al. 2006). 
The survey was conducted from two 1.3-meter telescopes, 
located at Mount Hopkins, Arizona and Cerro Tololo, 
Chile, in order to provide coverage for both the north¬ 
ern and southern skies, respectively. The focal plane of 
each telescope was equipped with three 256 x 256 pixel 
arrays, with a pixel scale of 2" x 2". Each field on the 
sky was covered six times, with dual 51-millisecond and 
1.3-second exposures, achieving a 10 it point-source depth 
of approximately 15.8, 15.1 and 14.3 mag (Vega) in J, H 
and K s , respectively. Calibration of the survey is con¬ 
sidered accurate at the 0.02 mag level, with photometric 
uncertainties for bright sources below 0.03 mag (Skrut- 
skie et al. 2006). 

The addition of 2MASS data requires an expansion 
of the stellar model described in Green et al. (2014) to 
cover the 2MASS .7, II and K s passbands, as well as an 
estimate of the survey selection function for 2MASS. We 
address these additions to our model in Appendices A 
and B. 

3.3. Source Selection 

We match each PS1 source to the nearest source in the 
2MASS point-source catalog, rejecting stars at greater 
than 2" separation. We require detection in at least four 
passbands, two of which must be PSl passbands. In 
order to reject extended sources, we require that m ps f — 
'"aperture < 0.1 mag in at least two PSl passbands. We 


TABLE 1 
PlXELIZATION 


nside 

Max. 

Stars / Pixel 

Solid Angle 
at Resolution (deg 2 ) 

# of 
Pixels 

64 

200 

77 

91 

128 

250 

90 

430 

256 

300 

11980 

228373 

512 

800 

16071 

1225471 

1024 

1200 

2957 

901971 

2048 

— 

66 

80956 

total 

— 

31240 

2437292 


additionally reject sources flagged as extended sources in 
2MASS. 

Individual PSl passbands are rejected if they have 
photometric uncertainty greater than 0.2 mag, or if the 
sources are beyond or close to the saturation limit for 
PSl (here considered 14.5, 14.5, 14.5, 14 and 13 mag 
in grizypi , respectively). In 2MASS passbands, we 
make the recommended “high-reliability catalog” selec¬ 
tion cuts, 8 in addition to the following requirements: 

• contamination/confusion flag = 0, 

• galaxy contamination flag = 0. 

Our final catalog contains 798,611,689 sources, of 
which 32% are detected in four passbands, 49% in five 
passbands, and 19% in six or more passbands. 

3.4. Pixelization 

We divide the sky into IIEALPix pixels (Gorski et al. 
2005), adjusting the pixel scale in order to keep the num¬ 
ber of stars per pixel roughly constant (see Fig. 2). Our 
procedure is to begin with nside = 64 pixels, and then 
subdivide each pixel recursively, as long as the number 
of stars exceeds some threshold, dependent on the pixel 
scale. We use thresholds given in Table 1, chosen to allow 
us to reach a resolution of nside = 512 with a relatively 
small number of stars, but to avoid going to higher res¬ 
olutions unless the stellar density is much higher. We 
reject pixels with fewer than ten stars. Such pixels com¬ 
prise a negligible fraction of the sky. In all, we assign just 
under 800 million stars, covering just over three quarters 
of the sky, to 2.4 million pixels, with an average of 327 
stars per pixel. 

4. RESULTS 

Applying the method described in Green et al. (2014), 
with the modifications described above, to Pan-STARRS 
1 and 2MASS stellar photometry, we produce a three- 
dimensional map of dust reddening, covering the three- 
quarters of the sky north of 8 = —30°. 

4.1. Distance Slices of Map 

In Figs. 4 and 5, we present the differential reddening 
in four spherical shells, centered on the Sun. Each panel 
shows the median dust reddening in a different range of 
Solar-centric distances. Due to perspective, dust at high 

8 See the 2MASS All-Sky Data Release Explanatory Supple¬ 
ment: http:// mww. ipac.caltech.edu/2mass/releases/allsky/ 

doc/secl_6b.html#composite 
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Fig. 4.- Median differential reddening in two Solar-centric distance ranges. The distance breaks are chosen to coincide with distance 
moduli // = 7.5 and 9.5, which line up with edges of distance bins in our map. We adopt a square-root stretch, in order to capture both 
low- and high-reddening features. The hole at £ as 120°, b as 30° corresponds to declinations above ~84°, which had not yet been fully 
processed at the time we created our 3D map. 


Galactic latitudes resides nearby, as Galactic dust lies in 
a thin disk. We recover the wealth of struct ure seen in the 
ISM across a wide range of scales, from thin filaments to 
large cloud complexes. Large, coherent cloud complexes 
appear at consistent distances. 

Fig. 6 gives a closer view of the anticentral region 
(l ~ 180°). Different features appear clearly in each ren¬ 
dered distance slice. The Perseus, Taurus and Auriga 
cloud complexes dominate the anticentral region in the 
closest distance slice, while the Orion molecular cloud 
complex (('. ~ 210°, b ~ —15°) and the California nebula 
(('. ~ 160°, b ~ —8°) appear very strikingly in the second 
distance slice. Of particular interest is the ring-like struc¬ 
ture that Orion A and B appear to be embedded within. 
This ring-like structure is only apparent when the back¬ 
ground dust is removed. In particular, the northeast por¬ 
tion of the ring is confused with the plane of the Galaxy 
in projection. Schlafly et al. (2014a) discusses the “Orion 


ring,” including possible formation scenarios for the ring, 
in greater depth. 

Fig. 7 shows the Galactic plane from l = 60° to 155° in 
more detail. The Cepheus flare, which lies at the center 
of the frame, at 95° < ( < 110°, b = 15°, separates 
into two components at different distances, visible in the 
first and third panels. Using a modified version of the 
method used in this paper along individual sightlines, 
Schlafly et al. (2014b) places the two components of the 
Cepheus flare at distances of 360 ± 35 pc and 900 ± 90 pc. 
The dust associated with Cygnus X (i % 80°, b as 0°) 
appears in the third panel, along with a wealth of fine 
structure along the Galactic plane. 

We note that each pixel is fit independently, and our 
only prior assumption about the spatial structure of 
the dust is that if forms a diffuse disk, as shown in 
Fig. 1. The detailed spatial structure in the interstel¬ 
lar medium that our analysis reveals indicates that the 
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Fig. 5.- As in Fig. 4, but with breaks at distance moduli n = 11.5 and 14.5. 


PSl and 2MASS photometry dominates over our priors 
out to a distance of several kiloparsecs and reddening 
of E (B — V) ss 1.5mag. With the assumption of spatial 
correlations between neighboring pixels, we expect that 
one would be able to significantly reduce the uncertainty 
in the map, and achieve better distance resolution (see, 
e.g., Sale & Magorrian 2014). 

4.2. Maximum & Minimum Reliable Distances in Map 

Our 3D dust map is based on measurements of stel¬ 
lar distances and reddenings. Beyond the most distant 
stars, we have no sensitivity to dust, and in front of the 
nearest stars, we have no information about the distance 
to the dust. Therefore, we estimate the minimum and 
maximum distance to which our map is reliable by locat¬ 
ing the nearest and farthest stars in each pixel. Outside 
of this distance range, our line-of-sight reddening infer¬ 
ences are dominated by our priors. Using the improved 
stellar parameter inferences (described in §2.4), we de¬ 
fine the minimum reliable distance in each line of sight 
as the distance out to v r hich there are iV c i oser observed 


main-sequence stars, and the maximum reliable distance 
as the distance beyond which there are A r f art her observed 
main-sequence stars. 

For this calculation, we count each Markov Chain sam¬ 
ple in stellar distance as a fraction of a star. We exclude 
stars that fail to converge, for which the model is a very 
poor match to the data (as determined by the Bayesian 
evidence for the star; see Green et al. 2014), or which are 
inconsistent with the inferred line-of-sight reddening pro¬ 
file. We consider a star inconsistent with the line-of-sight 
reddening inference if none of the 100 stored Markov 
Chain samples of the stellar distance and reddening is 
within a fractional distance crs (the modeled intra-pixel 
scatter in the dust column; see §2.4) of the line-of-sight 
reddening profile. Such objects are likely not well fit by 
any of our stellar templates, or alternatively signal that 
there is more variation in the dust column at fine angular 
scales than we allow. 

In determining the minimum and maximum reliable 
distances in the 3D dust map, we use only main- 
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Fig. G.- A closer view of the dust in the anticentral region in three distance bins. As in Figs. 4 and 5, we plot the median differential 
reddening. The Taurus-Perseus-Auriga complex is visible in the right half of the nearest distance bin. In the second distance bin, the Orion 
complex is visible on the left, while the California cloud is visible on the right. Note the ring-like shape of the Orion complex, which is 
only revealed by 3D mapping when confusion from background dust is removed. See Schlafly et al. (2014a) for discussion of this feature. 
Monoceros R2 appears beyond Orion, in the third distance bin, flanked by the plane of the Galaxy at yet greater distance. 
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sequence stars. This is because we consider our infer¬ 
ences for giants to be less reliable than our inferences 
for dwarfs. The colors and luminosities of giants depend 
more strongly on metallicity and age, the latter of which 
we do not model. 

In the left two panels of Fig. 8, we show the results 
obtained by requiring N c \ oser = 2 and iVf ai -ther = 10. The 
results are qualitatively similar for other choices of these 
parameters, as long as they are small compared to the 
typical number of stars in a pixel. The closest reliable 
distance is set almost entirely by the angular pixel scale. 
The nearby density of stars is, to a very rough estimation, 
uniform, meaning that the distance to the closest star in 
a pixel is a function primarily of the solid angle covered 
by the pixel. Accordingly, the top left panel of Fig. 8 is 
essentially a map of pixel solid angle, with boundaries in 
distance following pixel scale boundaries. 

The farthest reliable distance of the 3D dust map is 
strongly influenced not only by the pixel scale, but also 
by the distribution of stars throughout the Galaxy, the 
3D distribution of dust and the survey depth. The survey 
depth plays a larger role in the far limit because the most 
distant observed stars lie preferentially near the limiting 
survey magnitude. The closest stars, by contrast, are 
distributed more evenly in apparent magnitude. Accord¬ 
ingly, the high-latitucle patches with particularly shallow 
maximum reliable distances (e.g., the shallow patch near 
t = 205°, b = —55°) correspond to areas which are miss¬ 
ing one or more PS1 bands, or which have fewer obser¬ 
vation epochs in the PS1 3 tt survey, and therefore worse 
coverage. 

As a check on these estimated nearest and farthest re¬ 
liable distance maps, we construct equivalent maps using 
simulated stellar catalogs. For each pixel in our map, we 
use our Galactic priors and the estimated limiting mag¬ 
nitudes for grizypi and 2MASS JHK S to generate an 
equal number of stars as actually observed. We then de¬ 
termine the distance in front of which there are N c i osei 
and beyond which there are A’f art her main-sequence stars. 
In order to remove the complicating factor of the line-of- 
sight dust inference, we restrict this test to \b\ > 15° and 
to pixels for which E(B — F) r . r , < 0.05 mag, and assume 
for our simulated catalogs that there is zero dust extinc¬ 
tion. In these comparisons, we find that a halo number 
density close to the value given in , Iuric et al. (2008) is 
required to reproduce inferred map depths of Fig. 8. The 
results for a halo strength of n/, = 0.004, binned down to 
nside = 64, are shown in the right two panels of Fig. 8. 

These results indicate that while we dialed down the 
halo strength in our priors, the data, as reflected in 
our photometric stellar inferences, nonetheless prefers a 
stronger halo. In future work, we will investigate the im¬ 
plications of our stellar and line-of-sight dust inferences 
for a global Galactic structure model. 

4.3. Stellar Inferences 

As described in §2.4, we reweight, the naively inferred 
parameters for each star in order to take the line-of- 
sight reddening profile into account. In order to demon¬ 
strate the improvement in stellar inferences after this re- 
weigthing procedure, we compare our reddening infer¬ 
ences to a set of independently determined stellar red¬ 
dening standards, as in Green et al. (2014). For this, 
we use the set of stellar reddening standards developed 


by Schlafly & Finkbeiner (2011). The Sloan Extension 
for Galactic Understanding and Exploration (SEGUE; 
Yanny et al. 2009), part of SDSS-II, used moderate- 
resolution spectra to classify 240,000 stars. Empirically 
adjusted SDSS colors based on model atmospheres can 
then be compared with observed SDSS photometry to 
obtain reliable reddening estimates. We select the same 
sample of 200,000 SEGUE target stars as Schlafly & 
Finkbeiner (2011), which excludes white dwarfs and M 
dwarfs (the former because they are not contained in our 
model, and the latter because of their unreliable spectral 
classification). 

In Fig. 9, we compare reddening samples drawn 
from our PSl+2MASS-based stellar inferences to sam¬ 
ples drawn from the SEGUE-based estimates, which 
have Gaussian uncertainties. The agreement between 
the two reddening estimates improves significantly after 
reweighting our PS 1+2MASS-based stellar inferences. 
The improvement is most pronounced at low reddenings. 
The improvement is negligible for E {B — V) > 1 mag, 
due to the fact that we allow the individual stars to 
deviate from the local line-of-sight reddening profile by 
about 10%. Unlike the SEGUE-based reddenings, the 
PSl+2MASS-based reddening inferences are constrained 
to be non-negative, leading to a negative slope in the 
residuals near zero reddening. 

5. COMPARISON WITH PREVIOUS DUST MAPS 

Schlafly et al. (2014c) compares an earlier version of 
our 3D dust map with the two-dimensional, far-infrared 
emission-based SFD (Sclilegel et al. 1998) and Planck 
dust maps (Planck Collaboration et al. 2014). Here, we 
repeat this comparison using our new 3D dust map, and 
additionally compare our dust map with previous 3D 
dust maps, which are also based on stellar photometry. 

5.1. 2D Dust Maps 

It is possible to obtain a 2D dust map from our 3D map 
by projecting out distance, i.e., by taking the cumulative 
reddening out to some large distance. Such a map is 
of particular interest for extragalactic astronomy, where 
Milky Way dust is essentially a foreground screen to be 
removed. Several all-sky 2D maps of dust reddening al¬ 
ready exist, among them Burstein & Heiles (1982), based 
on Hi emission and galaxy number counts, Sclilegel et al. 
(SFD; 1998) and two recent reddening maps derived from 
Planck data (Planck Collaboration et al. 2014), which 
model dust optical depth and temperature from far in¬ 
frared emission, and then calibrate a dust optical depth 
to reddening relation. Near-infrared stellar colors have 
also been used to determine dust reddening more directly, 
as in the NICE/NICER/NICEST family of algorithms 
(Lada et al. 1994; Lombardi & Alves 2001; Lombardi 
2009, respectively), and Rowles & Froebrich (2009). 

Schlafly et al. (2014c) compares a 2D projection of an 
earlier version of our 3D dust map with the widely used 
SFD reddening map, as well as the newer Planck redden¬ 
ing maps. We repeat a number of the same tests for our 
new 3D dust. map. 

We begin, however, with a comparison between the 
map presented in this paper and the map presented in 
Schlafly et al. (2014c). The most important differences 
between the 3D maps used here and in Schlafly et al. 
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Fig. 7.- A closer view of the dust in the direction of Cepheus and Polaris flares and the eastern portion of the Great Rift, including 
Cygnus X. The Cepheus flare (95° < i < 110°, b as 15°) separates into two components in distance, visible in the first and third panels. 
The dust associated with the Cygnus X region {(. as 80°, b as 0°) appears in the third panel. 
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Fig. 8.- The minimum and maximum reliable distances in the 3D dust map. The left two panels are based on the inferred distribution 
of stars along each line of sight, while the right two panels are based on simulated stellar catalogs generated from our Galactic priors. We 
consider the map reliable if there are at least two stars closer than the given distance, and at least ten stars beyond the distance. Pixel 
size, survey depth and completeness, stellar density throughout the Galaxy and the presence of dust all affect the distribution of observed 
stars along the line of sight, and therefore the distance range over which the map is reliable. Sharp transitions in reliable distance occur 
at pixel size boundaries, where the number of stars per pixel changes discontinuously, as well as in patches of the sky with missing PS1 
passbands. 
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Fig. 9.— Histograms of the residuals of our PSl+2MASS-based 
stellar reddenings (referred to as “GSF” here), versus reddening 
estimates obtained by comparing SEGUE spectral classifications 
with SDSS photometry. The histograms are spread out along the x- 
axis by the local SFD reddening. The blue lines trace the 16th, 50th 
and 84th percentiles of the residuals in each bin of SFD reddening. 
In the upper panel, we use unweighted samples drawn from the 
individual stellar parameter Markov Chains. In the lower panel, 
we use reweighted samples, as described in the text. 


(2014c) are addition of near-infrared 2MASS photometry 


in our newer map, and that we sample here from the full 
posterior distribution on line-of-sight reddening, rather 
than finding the maximum-likelihood line-of-sight red¬ 
dening. Because Schlafly et al. (2014c) requires that each 
star be detected in gp i, and because we incorporate near- 
infrared 2MASS photometry alongside PSl photometry, 
our map reaches to deeper dust extinctions. This is ap¬ 
parent in the right panel of Fig. 10, where our new map 
predicts more reddening both in the inner Galactic plane 
and in dense dust clouds off the plane, such as Orion, 
Taurus and Perseus. 

At reddenings below E (B — V) < 0.08 mag, our new 
map predicts significantly less reddening than Schlafly 
et al. (2014c), as can be seen in the left panel of Fig. 
10. Beyond E (B — V) ~ 0.1 mag, our reddening scales 
agree very closely, with our new map predicting ~ 4% 
more reddening. The scatter between the two reddening 
measures comes to ~ 0.04 mag, with a maximum median 
offset of 0.04 mag at E(B-F) a; 0.08 mag, decreasing 
gradually to an offset of less than 0.01 mag between the 
two measures at E(J3 — V) = lmag. The behavior of 
the residuals suggests that at very low reddenings, our 
new 3D dust map prefers essentially zero reddenings too 
strongly. This may be due to the stronger priors on dust 
reddening used in the present work, or due to slight dif¬ 
ferences in onr new combined PS1-2MASS stellar locus, 
versus the PSl stellar locus used in Schlafly et al. (2014c). 

Next, we compare our inferred cumulative reddening 
out to 5 kpc with the SFD map, the Planck 7353 ghz- 
based reddening (hereafter, denoted simply as T 353 , 
with units of magnitudes of E (B — V)), and the Planck 
radiance-based map (hereafter denoted TZ, likewise with 
units of magnitudes of E(B — V) ). All three of these maps 
are derived by modeling dust emission, and should there¬ 
fore have different types of systematic errors than our 

















14 



Fig. 10. A comparison of our new 3D dust map (“GSF”), integrated to 4.5 kpc, with the 2D map presented in Schlafly et al. (2014c). 
The left panel shows the difference between our map and Schlafly et al. (2014c) as a function of a third reddening measure with uncorrelated 
errors, the Planck T 353 GHz-based dust reddening map. The right panel maps the median residuals between our map and Schlafly et al. 
(2014c) across the sky. All units are in magnitudes of E(B — V). 





Fig. 11. A comparison of the 2D SFD reddening map with our 
PSl-based map (“GSF”), integrated out to 5 kpc. The top panel 
shows the SFD reddening over the footprint of the PS1 survey, 
clipped to 0.25 mag, while the bottom panel shows the residuals 
after subtracting off the median cumulative reddening out to 5 
kpc predicted by our 3D dust map. Both panels use the same 
color scale. 


stellar reddening-based dust map. 

The upper panel of Fig. 11 shows the SFD reddening 
across the footprint of our map, while the lower panel 
shows the residuals after subtracting off our 5 kpc map. 
Off the Galactic plane, the residuals are small, while 
larger residuals are found in the plane of the Galaxy, 
where there is significant dust past 5 kpc, and where 
PSl does not necessarily detect stars beyond all of the 
dust. Near the Galactic center, one noticeable feature is 
a blue halo, where we infer more dust than SFD. The 
residuals are correlated with dust reddening in this area, 
indicating a scale offset between the two maps, rather 
than a constant offset. 

The “blue halo” is again apparent if we compare our 
map to the the Planck T353 map. In Fig. 12, we compare 
our inferred reddening at 5 kpc with SFD, T 353 , and 1Z. 
The left panels map the residuals across the PSl survey 
footprint on a square-root stretch, emphasizing small- 
amplitude differences. In the left panels, the Planck- 
based maps have been scaled by a constant factor to 


match our inferred PSl reddening for |6| > 20°. 

As the “blue halo” occurs in the direction of the nearby 
Aquila Rift, it is possible that the cloud has anomalous 
dust properties. For example, the grain size distribution 
or composition may vary, so that Ry = 3.1 is not a good 
assumption in the Aquila Rift. Another possibile expla¬ 
nation for the “blue halo” is that our stellar inferences 
may be systematically biased towards greater reddenings 
in this direction due to limitations in our Galactic model. 
Our priors do not, for example, include a radial metallic- 
ity gradient in the disk components of the Galaxy, which 
could lead to biased metallicity estimates for stars to¬ 
wards the Galactic center. We leave this question for 
future investigation. 

In the right panels of Fig. 12, we plot the residuals of 
our inferred reddening to 5 kpc with SFD and the Planck 
emission-based maps as a function of reddening. Here, 
we do not scale the Planck maps by any factor to bring 
them into alignment with our map. To conduct this com¬ 
parison, we compare the emission-based maps to multiple 
random realizations drawn from the posterior on redden¬ 
ing to 5 kpc from our 3D dust map. We restrict our 
comparison to high-Galactic-latitude regions (|5| > 30°), 
and cut out the ecliptic plane (1 3 | < 20°), where im¬ 
perfectly subtracted Zodiacal light might contaminate 
the emission-based reddening maps. When comparing 
our PSl-based reddening with SFD, we place the Planck 
7353 GHz-based reddening on the ic-axis, as its errors are 
uncorrelated with both maps on the y- axis. When dis¬ 
playing the residuals between our PSl-based reddening 
and the Planck reddening maps, we place SFD along the 
x -axis for the same reason. 

For reddenings above E(/J — V) > 0.05 mag, we see 
broadly similar residuals as found in Schlafly et al. 
(2014c), with with different behaviors above and below 
E(B — V) r: 0.15mag. Above E (B — V) « 0.15 mag, our 
map predicts about 10% less reddening than SFD, but 
is in good agreement with T 353 . As in Schlafly et al. 
(2014c), we And an overall difference in scale between 
our PSl-based map and the Planck 7?.-based map. For 
E(B — V) < 0.05 mag, we see the same residual between 
our map and the emission-based maps as found between 
our map and Schlafly et al. (2014c). For these small red¬ 
denings, our map favors essentially zero reddening too 
heavily. 
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Fig. 12. A comparison of our PSl-based reddening to 5 kpc (denoted by “GSF”) with SFD, the Planck T 353 GHz-based reddening, and 
the Planck radiance-based map (denoted by 1Z). The left panels map the residuals (between each map and the median integrated reddening 
in our map) on a square-root stretch. The panels on the right show the histogram of the residuals (between each emission-based map 
and posterior samples from our map) as a function of E (B — V). In the upper right panel, we put the T 353 GHz-based reddening estimate, 
in magnitudes, along the x-axis, since its errors should be largely uncorrelated with the errors in the PS1 and SFD reddening estimates. 
In the bottom two panels on the right, we use the SFD reddening estimate, in units of magnitudes, as our proxy for reddening, likewise 
because its errors should be uncorrelated with those of the quantities along the iy-axis. An inset in the top-right panel shows the regions 
that are masked in this analysis. The detailed behavior of the residuals, particularly at large reddenings, depends on which regions are 
masked, indicating that there are systematic differences in the residuals between our reddening map and emission-based map in different 
regions of the sky. 



The exact behavior of the residuals depends on which 
regions of the sky are masked in the analysis, indicat¬ 
ing that there are spatiallv-dependent systematic differ¬ 
ences in the residuals between our reddening map and 
emission-based maps. However, the essential features 
remain the same, with different slopes in the residuals 
below and above E (D — V) ~ 0.15mag, and our map fa¬ 
voring lower reddening below E (B — V ) ~ 0.05 mag. 

5.2. Marshall et al. (2006) 

Marshall et al. (2006) developed a method to determine 
the reddening-distance relation along individual lines of 
sight by comparing 2MASS J — K s stellar colors to those 
of simulated catalogs based on the Besangon model of the 
Galaxy (Robin et al. 2003). Marshall et al. (2006) then 
applied this method to a regular grid of sightlines sepa¬ 
rated, by 15' covering the region |£| < 100°, \b\ < 10°. 
The result is a 3D map of reddening in the inner Galaxy, 
extending to a maximum extinction of ~ 1.4 — 3.75 mag 
in the 2MASS K s band (equivalent to ~ 4.5 — 12 mag in 
E(B — V')), and a maximum distance of ~ 7kpc. Because 
Marshall et al. (2006) use only giants in their analysis, 
the dust map they produce has little information in the 
nearest kiloparsec. We will refer to this map as the “Mar¬ 
shall map.” 

Our dust map overlaps with the Marshall map in 
the approximate region 0° < l < 100°, |6| < 10°. 


Fig. 13 shows the cumulative reddening at increasing 
distances in both the Marshall map mid our 3D dust 
map. When converting from extinction to reddening, 
we assume A Ks = 0.320 x E (B — V), as calculated by 
Yuan et al. (2013) for a 7000 K source spectrum at 
E(B-V) = 0.4mag, using the Cardelli et al. (1989) red¬ 
dening law and assuming By = 3.1. We mask our map 
beyond our predicted maximum reliable distance, as de¬ 
termined in §4.2. In the large-distance limit, our maps 
show good qualitative agreement outside of the masked 
areas. Fig. 14 shows the differential reddening in the two 
maps in bins of increasing distance. 

In the nearest two to three kiloparsecs, our map shows 
clearly differentiated structures at discrete distances that 
are spread over several distance bins in the Marshall dust 
map. For example, the Cygnus rift (located at i ~ 80°, 
b ~ 0 °) appears clearly in our map in the distance bin 
spanning 0.5 — 1 kpc, while in the Marshall map, it is 
spread over all the distance bins closer than ~ 3 kpc, 
and is therefore not cleanly separated from superimposed 
dust structures at greater distances. The greater distance 
resolution of our map at nearby distances is most likely 
due to the fact that we use both main-sequence stars and 
giants, while Marshall et al. (2006) relies solely on giants, 
which are saturated nearby and form a larger fraction of 
the observable stellar population at greater distances. 

In addition to better distance resolution in the first 
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Fig. 13. A comparison of the median cumulative reddening out to increasing distances in the Marshall map (left panel) and our map 
(right panel). Regions beyond the maximum reliable distance in our map are masked out in blue in the right panel. At large distances, 
the two maps agree qualitatively, with the masked regions in our map corresponding to the most heavily obscured regions in the Marshall 
et al. (2006) map. The Marshall map has greater depth, but lower angular resolution, and lower distance resolution in the nearest two to 
three kiloparsecs. 


few kiloparsecs, the greater source density of PS1 rela¬ 
tive to 2MASS allows us to achieve better angular reso¬ 
lution than the Marshall map. Fig. 15 demonstrates the 
difference in angular resolution between the two maps. 
In most of the region of overlap between our reddening 
map and the Marshall map, we achieve an angular reso¬ 
lution of 3.4', as opposed to the constant 15' resolution 
of the Marshall map. This allows us to resolve detailed 
filamentary structure not seen in the latter map. As 
dust reddening can vary significantly on small angular 
scales, this increased angular resolution will be impor¬ 
tant in practice for correctly de-reddening extinguished 
sources in regions with complex dust structure. 

Deep in the plane of the Galaxy, where high extinc¬ 
tion reduces source counts, our finer angular resolution 
limits the distance to which our map can trace dust, com¬ 
pared to the Marshall map. Although the PS1 3 tt sur¬ 
vey is deeper than 2MASS, the 2MASS passbands are 
less affected by dust extinction, and the advantage of 
PSl decreases in regions of high extinction. In such re¬ 
gions, the greater depth of PSl does not fully compen¬ 
sate for our smaller pixels, limiting the depth to which 


we trace dust deep in the Galactic plane. The fact that 
our map derives its most accurate reddening informa¬ 
tion from main-sequence stars also limits the maximum 
extinction to which it is reliable. 

5.3. Lallement et al. (2013) 

Combining distance and reddening estimates for ~ 
23,000 stars with the assumption of spatial correlation 
in dust density, Lallement et al. (2013) infer reddening in 
3D out to a distance of ~ 800 pc from the Sun. We find 
close morphological agreement between our 3D dust map 
and that of Lallement et al. (2013), with some differences 
which are worth taking note of. 

In Fig. 16, we show the distribution of dust and stars in 
a slice 25 pc above the Galactic plane, level with the Sun. 
We show the median dust, density. In order to generate 
the stellar locations, we draw a sample at random from 
the improved distance posterior of each star (see §2.4), 
and select only those stars that lie within 5pc of the 
chosen plane. For display purposes, we only display one 
out of every thousand stars. 

Just as Lallement et al. (2013), we find cavities in the 
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Fig. 14.— A comparison of the median differential reddening in bins of increasing distance in the Marshall map (left panel) and our 
map (right panel). In the nearest two to three kiloparsecs, we achieve much better distance resolution, as evidenced by the differentiated 
structures visible in successive distance bins. 



Fig. 15.— A zoom-in of median cumulative reddening to 3kpc, showing the difference in angular resolution between the Marshall 3D 
reddening map (left panel) and our map (right panel), as well as the lower noise of the latter. The Marshall dust map has an anuglar 
resolution of 15', while our dust map has a typical angular resolution of 3.4' in the region of the Galactic plane displayed above. The 
reliability mask has not been applied our map in this figure. 
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dust density in the directions of £ = 70° and 225°. The 
overall morphology of the dust structure in the right 
planel of Fig. 16 matches that of Fig. 1 in Lallement 
et. al. (2013). 

The most obvious difference between our maps and 
those of Lallement et al. (2013) is the different voxel 
shapes employed in our work and theirs. Lallement 
et al. (2013) use small cubic voxels, and assume a spa¬ 
tial correlation function that favors similar dust densi¬ 
ties in nearby voxels. This allows them to densely sam¬ 
ple the reddening distribution, with voxels that are not 
directly constrained by stellar reddening measurements 
being constrained by neighboring voxels. In contrast, we 
infer the dust distribution in each line of sight separately, 
without assuming spatial correlations in the dust density, 
as laid out in Green et al. (2014). While our map has 
excellent angular resolution, it has distance bins with a 
width of about 25%, giving our voxels their pencil-beam 
shape. 

Although we see roughly the same structures, such as 
cavities in the reddening distribution centered on l = 70° 
and 225°, the distances we derive for a number of clouds 
is greater than the distances Lallement et al. (2013) find. 
In particular, while Lallement et al. (2013) place the 
Cygnus rift between 500 and 600 pc, we place it at a 
distance of 800 to 1000 pc. 

6 . ACCESSING THE MAP 

Our 3D dust map can be accessed at http:// 
argonaut. skymaps. inf o. The website provides an in¬ 
terface for querying individual lines of sight, as well 
as the ability to download the entire map and soft¬ 
ware to read it. We also provide an API through the 
website, which allows users to query the map remotely 
with a few lines of code, without the need to down¬ 
load the entire data cube. The data is also accessible 
at http://dx.doi.org/10.7910/DVN/40C44C, through 
the Harvard Dataverse. 

6.1. Data Cube 

The basic data product our map contains is samples of 
the differential reddening in 31 distance bins, in each of 
2.4 million pixels. The data structure is thus 

(2.4 x 10 6 pixels) x (500 samples) x (31 distance bins) . 

This allows us to determine the probability density of the 
cumulative reddening to any distance (within a few kilo- 
parsecs) in the 3 tt steradians covered by the PS1 survey. 

We encourage users to use the Markov chain samples 
of reddening versus distance provided by our interface in 
their analyses, as the samples contain the full statistical 
information generated by our method. These samples 
can be queried using our web API, downloaded as an 
ASCII table for individual lines of sight, or accessed di¬ 
rectly in the complete data cube provided for download. 

An additional, and larger, dataset that we produce 
while generating the 3D dust map is a library of photo¬ 
metrically determined stellar parameters. As described 
in Green et al. (2014), we infer the probability distri¬ 
bution of distance modulus, reddening, metallicity and 
absolute rpi-magnitude for each star. We thus have a 
second data cube with shape 

(8 x 10 8 stars) x (100 samples) x (4 parameters) . 


In addition, we store quality assurance information for 
each star, including whether the Markov Chain con¬ 
verged during the fitting procedure, and the Bayesian 
evidence for the stellar model, which is similar to the 
X 2 statistic in maximum-likelihood fitting. Point sources 
with poor evidence are likely either of stellar types not 
contained in our model, such as very young stars or 
binary systems, or are not stars (e.g., white dwarfs, 
quasars, unresolved galaxies). 

7. CONCLUSIONS 

We have presented a three-dimensional map of dust 
reddening covering three quarters of the sky, based on 
photometric inferences for ~ 800 million stars with high- 
quality multiband photometry in Pan-STARRS 1. The 
map provides a window into the structure of the in¬ 
terstellar medium, revealing detailed structure from the 
smallest scales in our map, 3.4', all the way to large 
cloud complexes spanning many degrees. We provide 
interfaces to query and download the map at http: 
//argonaut.skymaps.info. 

Projecting our map down to two dimensions, we find 
good agreement with emission-based dust maps at high 
Galactic latitudes, where we expect our method to trace 
the entire dust column. Comparison with the 3D map of 
nearby dust presented in Lallement et al. (2013) shows 
the same morphological features. In comparison with the 
3D dust map of the inner Galactic plane presented in 
Marshall et al. (2006), our map has greater angular res¬ 
olution and superior distance resolution within ~ 3 kpc. 
However, due to our reliance primarily on main-sequence 
stars, as opposed to giants, our map does not penetrate 
to as great a depth as Marshall et al. (2006). 

Our dust mapping technique can be extended to 
take advantage of photometric surveys beyond PS1 and 
2MASS. The Dark Energy Survey (DES; The Dark En¬ 
ergy Survey Collaboration 2005) is surveying 5000 deg 2 
of the southern sky, largely complementary to the PS1 
footprint, in a similar filter set. The LSST (Ivezic et al. 
2008) will provide deep ugrizy photometry for the sky 
south of S ~ 34.5°. In the nearer future, the ESA Gaia 
mission (Lindegren et al. 1994) will provide multiband 
photometry, geometric parallax distance measurements 
and proper motions for one billion stars. Parallax dis¬ 
tances, where available, will vastly improve stellar dis¬ 
tance estimates, breaking the dwarf-giant degeneracy in 
particular. Kinematic information from Gaia will pro¬ 
vide additional information about, stellar distances and 
types. 

A reliable map the 3D distribution of dust is impor¬ 
tant to studies of stellar populations within the Galaxy, 
as well as streams and global morphology of the Milk 
Way. In future work, we will leverage the stellar infer¬ 
ences produced as a by-product of our map to study the 
morphology of the Galaxy. We expect that the three- 
dimensional dust map presented here will find many dif¬ 
ferent uses not yet envisioned by the authors. 
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Fig. 16. Dust density in a lOpc-thick slice lying 25 pc above the Galactic plane (i.e., level with the Sun), with positions of stars within 
this slice overplotted in the left panel. In detail, we show the median (over multiple realizations of the 3D dust map) of the reddening 
column density along each sightline passing through the slice (i.e., into the page). The Sun is at the origin of each panel, with the right 
panel being a zoomed-in version of the left (without the overplotted stars). Only every 1000th star has been plotted. The dust map is 
reliable out to distances at which the stellar density becomes too small to trace the dust column. Compare with Fig. C2 of Lallement et al. 
(2013), in which the stellar density is concentrated within ~200 pc of the Sun. 


ipating institutes, the Max Planck Institute for Astron¬ 
omy, Heidelberg and the Max Planck Institute for Ex¬ 
traterrestrial Physics, Garching, The Johns Hopkins Uni¬ 
versity, Durham University, the University of Edinburgh, 
Queen’s University Belfast, the Harvard-Smithsonian 
Center for Astrophysics, the Las Cumbres Observatory 
Global Telescope Network Incorporated, the National 
Central University of Taiwan, the Space Telescope Sci¬ 
ence Institute, the National Aeronautics and Space Ad¬ 
ministration under Grant No. NNX08AR22G issued 
through the Planetary Science Division of the NASA Sci¬ 
ence Mission Directorate, the National Science Founda¬ 
tion under Grant No. AST-1238877, the University of 
Maryland, and Eotvos Lorand University (ELTE). 

This publication makes use of data products from the 
Two Micron All Sky Survey, which is a joint project of 
the University of Massachusetts and the Infrared Pro¬ 
cessing and Analysis Center/California Institute of Tech¬ 
nology, funded by the National Aeronautics and Space 
Administration and the National Science Foundation. 

The computations in this paper were run on the 
Odyssey cluster supported by the FAS Science Division 
Research Computing Group at Harvard University. 

Gregory Green and Douglas Finkbeiner are supported 
by NSF grant AST-1312891. Edward Schlafly acknowl¬ 
edges funding by Sonderforschungsbereich SFB 881 “The 
Milky Way System” (subproject A3) of the German Re¬ 
search Foundation (DFG). Nicolas Martin gratefully ac¬ 
knowledges the CNRS for support through PICS project 
PICS0G183. 

APPENDIX 

PAN-STARRS 1 AND 2MASS STELLAR 
TEMPLATES 

Green et al. (2014) describes how stellar templates are 
compiled for the PS1 passbands. Briefly, metallicity- 
independent main-sequence stellar colors were obtained 
by fitting a stellar locus in color-color space, and 


metallicity-dependent absolute magnitudes were ob¬ 
tained from the metallicity-dependent photometric par¬ 
allax relation given in Ivezic et al. (2008). For the giant 
branch, linear fits to globular cluster color-magnitude di¬ 
agrams from Ivezic et al. (2008) were used. 

In order to include 2MASS photometry in our dataset, 
we require joint PS1-2MASS stellar templates. We com¬ 
pile these templates using a. nearly identical procedure 
as in Green et al. (2014). We begin by selecting ~1 mil¬ 
lion stars with E(B— U) SFD < 0.1 mag, detections in all 
PSl and 2MASS passbands, and photometric errors less 
than 0.5 mag in every passband. The resulting sample 
has a median reddening of 0.016 mag in E (B — V). After 
dereddening the photometry, we fit a stellar locus in 7- 
dimensional color space, using the algorithm laid out in 
Newberg & Yanny (1997). The resulting stellar locus is 
plotted in Fig. 17. 

As before, we apply a metallicity-dependent photomet¬ 
ric parallax relation to obtain a set of stellar templates, 
indexed by absolute rpi-magnitude and [ Fe / H ]- In order 
to obtain templates for the giant branch which incor¬ 
porate 2MASS, we again use the templates from Ivezic 
et al. (2008), replacing stellar colors with our new joint 
PS1-2MASS stellar locus colors. We use rpj — i p i to 
match each giant template to a color template in our 
PS1-2MASS stellar locus. 

In the 2MASS passbands, we adopt the reddening coef¬ 
ficients R.j = 0.786, Bh = 0.508 and Bks = 0.320, calcu¬ 
lated by Yuan et al. (2013) for a 7000 K source spectrum 
at E(B — V) = 0.4 mag, using the Cardelli et al. (1989) 
reddening law and assuming By = 3.1. 

2MASS SELECTION FUNCTION 

In order to avoid Malmquist bias in photometric infer¬ 
ences, it is necessary to characterize the survey selection 
function. In Green et al. (2014), we characterized this as 
the detection probability in each passband as a function 
of apparent magnitude and position on the sky. Denote 
the detection or non-detection of a point source on a 
particular location on the sky in passband i as a binary 
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Fig. 17.— Stellar locus fit in join PS1-2MASS color space. For this fit, ~1 million star’s in the vicinity of the North Galactic Pole, with 
SFD reddening less than 0.1 mag in E(B— V), were used. The fitted stellar locus anchor points are plotted over the density of stars in each 
color-color projection. A reddening vector with magnitude 0.5 mag in E (B — V) is overplotted for each combination of colors. 


variable, S t . Denote the intrinsic (or modeled) appar¬ 
ent magnitude of the star in passband i as m mo d,i- The 
detection probability is written as 

p{Si | ?n m od.i) • (Bl) 

For point-source detections in PSl, we modeled the de¬ 
tection probability based on the sky background, read 
noise and point-spread function. As the 2MASS point- 
source catalogs do not contain this information, we model 
the 2MASS selection function empirically, from the his¬ 
togram of detections as a function of apparent magnitude 
across the sky. 

There are two ways to approach this problem. The first 
approach is to construct a full forward model, drawing 
stellar types, locations and reddenings from our Galactic 
and stellar priors, generating model photometry for the 
simulated stars, and applying a trial selection function 
and photometric errors to obtain a sample of simulated 
“observed” apparent magnitudes. One would then vary 


the selection function until the histogram of detections 
versus observed apparent magnitude matched observed 
histogram in a given region of the sky. This method 
suffers from its reliance on our relatively crude priors on 
Galactic reddening, and its sensitivity to any errors in 
our Galactic and stellar priors. 

We therefore opt for a simpler approach. We make 
the assumption that near the limiting magnitude, the 
true sky density of objects is a smooth function of ap¬ 
parent magnitude. In small patches of the sky, we lo¬ 
cate the turnoff in the histogram of detections, m toi , 
as a function of apparent magnitude. In the range 
—3.5 < nij — m to , * < —0.1, we model the logarithm of 
the number of detections in each apparent magnitude bin 
as a third-order polynomial in magnitude. This is our 
smooth estimate of the intrinsic sky density of sources 
as a function of apparent magnitude. For each 2MASS 
passband, we define the limiting magnitude as the bin in 
which the observed sky density per unit magnitude falls 
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Fig. 18. Limiting magnitude of the 2MASS survey in all three 
passbands, across the PS1 survey footprint. The limiting magni¬ 
tude is defined as the apparent magnitude at which 50% of point 
sources are detected. Effective survey depth is relatively uniform 
over most of the sky, but falls off precipit ously towards the Galactic 
Center due to source crowding. 


below 50% of the estimated true sky density per unit 
magnitude. 

We use this procedure to construct maps of the 2MASS 
J, II and K s limiting magnitudes at two HEALPix res¬ 
olutions, nside = 32 and 64. In each region of the sky, 
we adopt the nside = 64 map, unless the given pixel 
contains fewer than 1000 detections in the given pass- 
band, in which case we switch to the lower resolution, 
nside = 32 map of limiting magnitude. Fig. 18 shows 
the resulting, multi-resolution map of limiting magnitude 
in each passband. 

MARKOV CHAIN “SWAP” PROPOSALS 

Because of the way in which we parameterize the line- 
of-sight distance versus reddening profile, there can be 
strong anticorrelations between the different parameters. 
The differential reddening in one distance bin is often 
strongly anticorrelated with the differential reddening in 
neighboring distance bins, because transferring some dif¬ 
ferential reddening from one bin to a neighboring bin 
keeps the cumulative reddening constant in following dis¬ 
tance bins. Similarly, if the data is well fit by a large 
jump in reddening in one bin, it is often also well fit 
by a jump in reddening in a neighboring bin. These 
strong anticorrelations between the differential redden¬ 
ing in neighboring distance bins can slow down Markov 
chain convergence when fitting the line-of-sight redden¬ 
ing profile. 

We therefore introduce a new type of Markov Chain 
Monte Carlo (MCMC) proposal step, which we term the 
“swap” proposal. In order to generate a proposal state 



Fig. 19. Generation of an MCMC proposal state by the “swap” 
proposal method. The proposal state is identical to the current 
state, up to a swap in differential reddening in two bins. The black 
curve shows the cumulative reddening of the current state, while 
the light gray curve shows the cumulative reddening of the proposal 
state. The blue bars show the current differential reddening in 
each distance bin, while the striped green bars show the proposed 
differential reddening in each distance bin. Note that the two are 
identical, up to a swap between two distance bins. 


for the Maikov chain, we swap the differential reddening 
in two distance bins, as shown in Fig. 19. 

The swap proposal is symmetric, allowing the usual 
Metropolis-Hastings acceptance probability to be used. 
Call the current state A’ and the proposal state Y, and 
denote the two distance bins that were swapped to gen¬ 
erate Y as i and j. The probability of proposing Y, 
Q(X —> Y), is simply the probability that bins i and j 
are selected to be swapped. As long as the choice of i and 
j has nothing to do with the current state of the chain, 
the reverse step, from Y to A, is equally likely. That is, 
Q (Y —» A) is also simply equal to the probability that 
distance bins i and j are selected to be swapped. The 
acceptance probability, A(X —> Y), is therefore given by 


A{X -> Y) 


min ^1, 
min ^1, 


p(Y) Q(Y —> A) \ 
p(X) Q(X —» Y)J 
P(Y) \ 

P(X)J ’ 


(Cl) 


where p(X) and p(Y) are the probability densities of A 
and Y, respectively. This is simply the usual Metropolis- 
Hastings acceptance probability. 

The addition of swap proposals, alongside Metropolis- 
Hastings proposals and affine stretch proposals, allows 
the Markov chain to mix more quickly. Many distance¬ 
reddening curves for the sightline plotted in Fig. 3, for 
example, differ primarily in which distance bin the large 
jump in reddening occurs. In such a sightline, the addi¬ 
tion of swap proposals allows the Markov chain to transi¬ 
tion quickly between probable states, greatly improving 
convergence times. 


SCATTER IN THE LINE-OF-SIGHT 
REDDENING PROFILE 

As sketched out in 2.2, we allow each star to devi¬ 
ate slightly from the line-of-sight distance-reddening re¬ 
lation in the pixel. Allowing stars at the same distance 
to have somewhat different reddenings renders our model 
more robust to sub-pixel angular variations in dust den¬ 
sity. What follows below is a detailed derivation of how 
marginalization over S t . the deviation of star i from the 
average reddening profile in the pixel, translates into 
a. Gaussian smoothing of the stellar probability density 
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function in distance-reddening space, p(pi,Ei \ rhi). Be¬ 
ginning with the full posterior density on line-of-sight 
reddening given our stellar photometry, and including a 
deviation, 6 ,, for each star i, 

p{a | {m}) cx p(d) PJ I d pi d@i dii, p(?n, | p i: § { ,a, 

i 

x p(^p i: Qi,Si\a) . (Dl) 

Here, 0, is the type of star i. Taking just the integrand, 
and ignoring the subscript i for the moment, 

I = p(m | p, 0, a, <5^ p(^p, 0, 6 | oj (D2) 

= p(jh | p, 0, a, j) p(/i, 0 | a) | p, 0, cv) . (D3) 

As in Green et al. (2014), we assume that 

p(/r,0|a) =p(/r,0) , (D4) 

with the only complication being the survey complete¬ 
ness limit, which is dealt with explicitly in that paper, 
but which has no impact on the modification being dis¬ 
cussed here. We also assume that 

p(s | p, 0, o) = p(S | p, a) , (D5) 

since 6 is a property of the sub-pixel variation in dust 
density, and should be unrelated to stellar type. Then, 

I = p(m | p, 0, <5, p(p, ( r )) p{6 |p, a) . (DO) 


Recall that the prior on <5 is a Gaussian centered on zero, 
with width determined by E(p;3). Then, 

f dQ I x p(p,E = (1 + 6) E(p;a) \rh)p(S\E(p-,d)) , 

(Dll) 

Marginalizing over 6 now leaves us with a smoothed ver¬ 
sion of the individual stellar posterior density in distance 
and reddening: 

!dQd61 = p(p,E(p;a) \ rh) , (D12) 

where we have defined the “smoothed” individual stellar 
posterior probability density 

p(p. E) = J d<5 p(p, (1 -I- 6) E | rh) p(5 \ E) . (D13) 

In Eq. (3), we therefore replace the individual stellar 
posterior densities with “smoothed” posterior densities: 

p{a \ {rh,}) x p(a) J dpi p{p u E(pp a) \ rhi) . (D14) 

Our method relies on calculating the individual stellar 
posterior probability densities for all stars along a. line 
of sight first, before sampling from the line-of-sight red¬ 
dening distribution, as described in detail in Green et al. 
(2014). The above derivation shows that, this new addi¬ 
tion to our method, allowing stars to deviate from the 
line-of-sight reddening profile, requires an intermediate 
step, in which the individual stellar posterior densities 
are smoothed in reddening, according to Eq. (D13). 

STELLAR SAMPLE REWEIGHTING 


The likelihood term (the first term on the right-hand 
side), can be rewritten as 

p(m\p,Q,E(p-,a,5)') , (D7) 

since the modeled apparent magnitude is simply deter¬ 
mined by the stellar distance, type and reddening, and 
the individual stellar reddening is determined by the line- 
of-sight reddening profile, the stellar distance, and the 
fractional deviation of the stellar reddening from the lo¬ 
cal reddening. 

By Bayes’ Rule, the product 

p(rh\p,Q,E(p-,d,sjj p(p,©j (D 8 ) 

is proportional to the posterior density on distance, red¬ 
dening and stellar type for an individual star, in the pres¬ 
ence of a flat prior on reddening. Thus, 

I oap(^p,e,E\mj p(5\p,a) , (D9) 

evaluated at E = (1 + <5) E(p; a). After marginalizing 
over stellar type, 0 , we are left with 

j d© I x p(p, E) = (1 + 6) E(p\ a \ rh) p(6 \ p, a) , 

(DIO) 


Here, we derive in detail how to reweight the Markov 
Chain samples resulting from the naive stellar inferences, 
in order to condition the stellar parameters on the line- 
of-sight reddening profile. This discussion complements 

§ 2 - 4 ; 

Given a fixed line-of-sight reddening profile parameter¬ 
ized bv, E(p;a), each star is described by a distance p, 
stellar type 0, and fractional offset 6 from the reddening 
profile. The posterior of these parameters, conditioned 
on the star’s photometry and the line-of-sight reddening 
profile, is given by 

p(^p, 0 , 6 1 m, d^ oc p(jh \ p, ©, 5, a'j p(^p, 0 , 6 | oej . 

(El) 

The second term 011 the right-hand side can be broken 
down into two parts, 

p(/b $ I a) oc p{p, 0) P(S | p,a) . (E2) 

Thus, 

p(p, 0, 5 1 rh. x p(jh | p, Q. 6, dj p^p, Qj 

x p(S | p, a) . (E3) 

The two stellar parameterizations, |/r, 0, 6, cvj and 
| p. 0 . L'| , are equivalent, as the stellar distance mod- 
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ulus, p, the line-of-sight reddening parameters, a, and 
the star’s fractional offset, <5, from the line-of-sight red¬ 
dening are sufficient to determine the stellar reddening, 
E. We would like to reweight the samples we store in 

the space j p. 0, E j, so that they correspond to the pos¬ 
terior density given above, in Eq. (E3). As described in 
Green et al. (2014), the samples we store in our initial 
processing are drawn from the posterior density 

p(/b ©, E | m) oc p(m | p, 0, e) p(p, ©) p{E) , (E4) 


with a flat prior on E. so that p{E) = const. Trans¬ 
forming to the parameterization j/i, 0,(5, o|, we have 

to transform the flat prior p(E) to an equivalent prior, 
p(8 1 //., a). This prior is given by 


p{8 | p, 3) 


p(E) 


9_E_ 

98 


oc 


/l,Q 


9£ 

98 


li,a 


(E5) 


where we have taken out the constant factor p(E) on the 
r.h.s. and replaced the equality with a proportionality. 


The stellar reddening is related to p. a and 8 by Eq. (6). 
Taking the derivative w.r.t. 8, 

= Sc [(! + Si) E{pi ; a)} = E(pp a) . (EG) 

li,a 

The flat prior in E thus implies a prior on 8 equal to 

p{8 1 p. a) oc E(pi\a) . (E7) 

In order to sample from Eq. (E3), we can weight the 
samples from our initial processing by the ratio of the 
prior in the new parameterization, where reddening is 
conditioned on the line-of-sight reddening profile, to the 
prior used in the initial sampling, where a flat prior on 
reddening is used. From the above calculation, this ratio 
is given by 


9E_ 

98 


Pnew (8 Ip - a ) P new (£ | P , «) 

- (X - 

poid(<51 p, a) E{pi\a) 


(E8) 


The functional form of p new (<51 p, a) is given by Eq. (7). 
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